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Long wavelength oscillations (Tkachenko waves) of the triangular lattice of quantized vortices in 
superfluid neutron stars have been suggested as one of the possible explanations for the timing noise 
observed in many radio pulsars, in particular for the 100-1000 day variations in the spin of PSR 1828- 
11. Most studies to date have, however, been based on the hydrodynamics developed for superfluid 
Helium. In this paper we extend the formulation to a two fluid neutron and proton system, relevant 
for neutron star interiors and include the effect of chemical coupling, compressibility and mutual 
friction between the components. In particular we find that chemical coupling and compressibility 
can have a drastic effect on the mode structure. However, for the slower pulsars rotating at 1-10 
Hz (such as PSR B1828-11), most choices of parameters in the equation of state lead to Tkachenko 
oscillations with frequencies in the correct range to explain the timing noise. We also investigate the 
case of more rapidly rotating pulsars (above 100 Hz) for which we find that there is a vast portion 
of parameter space in which there are no Tkachenko modes, but only modified sound waves at much 
higher frequencies. 

I. INTRODUCTION 

A growing number of radio pulsars have now been observed for over a decade (some for more than 30 years) and are, 
in general, very stable rotators. However, many pulsars also exhibit timing irregularities, such as "glitches" , which 
are sudden increases in the rotation rate, and "timing noise" , a general term which refers to low frequency quasi- 
periodic structures that appear in the timing residuals, once the "regularly" pulsating solution has been removed. In 
particular, while the irregularities in younger pulsars are dominated by the recovery from glitch events, for a handful 
of older pulsars there is growing evidence for long-period (~ 100-1000 days) oscillations in the timing residuals [I]. 
In some cases these periodicities, and the correlated pulse shape changes, can be partially explained by neutron star 
free-precession, and one of the best examples of this is PSR B1828-11 which shows significant periodicity at ~ 256 
days and ~ 511 days [2J. However there are theoretical arguments stating that mutual friction between the interior 
superfluid components of the star would damp out any precessional motion on a short timescale [21 H] (although 
Glampedakis et al. [5] have shown that short wave length instabilities in the pinned superfluid could cast an element 
of doubt on such conclusions) . Furthermore recent work shows that some pulsars may be switching abruptly between 
two different states with different spin-down rates, thus giving rise to the observed timing behaviour [5J. 

Noronha and Sedrakian [7 , following earlier suggestions by Ruderman [5] , indicated that an alternative explanation 
for the observed long term periodicity could be the propagation of Tkachenko waves in the star. In fact, it has been 
suggested that Tkachenko waves excited by glitches may be driving precession in one of the X-ray Dim Isolated 
Neutron stars (XDINs), RX J0720.4-3125 0. 

Neutron star interiors are expected to contain charge neutral superfluids that rotate by forming an array of quantized 
vortices. In their lowest energy state the vortices form a two-dimensional triangular lattice that can support elastic 
oscillations, Tkachenko waves [10] . that have been studied extensively, both theoretically and experimentally, in 
superfluid 4 He (see e.g. [TT] for a review) and recently in Bose-Einstein condensates (BECs) (see e.g. [T2] and [T3]). 
The undamped propagation of Tkachenko waves in a neutron star would lead to periodic variations in the angular 
momentum of the superfluid which, due to coupling to the crust, would lead to variations in the observed rotation rate. 
In order to ascertain if this is a viable hypothesis it is crucial to understand how the detailed microphysics of neutron 
star interiors affects the propagation of the modes. Most studies to date have been based on the hydrodynamical 
theory of Tkachenko waves developed by Baym and Chandler [TH \TE\ for superfluid 4 He. In this case the fluids can 
be treated, to a good degree of approximation, as incompressible, given that the rotation rate is always well below the 
sound wave frequency (note, however, that in BECs, compressibility has a strong effect, due to the interactions being 
much weaker and the sound speed much lower than in Helium) , and the system can be described as a condensate (the 
"superfluid") coupled to a "normal" fluid which consists, loosely speaking, of the thermal excitations of the system. 

In a realistic neutron star, on the other hand, one must take into account not only the effects of rapid rotation, but 
also the presence of several massive fluids, describing the flow of electrons, protons, superfluid neutrons (and their 
excitations at finite temperature) and possibly exotic particles such as hyperons or deconfined quarks, which cannot, 
in general, be assumed to be incompressible. Furthermore, one must consider various dissipative processes that damp 
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out the oscillations. It is well known that in a multifluid system there will, in general, be many more dissipation 
channels than in a simple one- fluid flow described by the Navier-Stokes equation [16j E! • 

Solving the full problem is clearly a daunting task, so in this paper we shall make a series of simplifying assumptions. 
We take a system of two fluids: the superfiuid neutrons and a charged component of protons and electrons, which 
are locked by the Coulomb interaction on a much shorter timescale than the dynamical timescales considered here. 
This assumption is justified as the frequencies of the modes we calculate (Tkachenko waves, sound waves and inertial 
waves) are always significantly lower than the electron-proton cyclotron and plasma frequencies, i.e. lower than w 10 15 
Hz (|18|). Both fluids are assumed to be compressible, and we shall assume some simplified analytic models for the 
equation of state, derived from Haskell et al. |19j . Furthermore we will only consider the damping due to superfiuid 
mutual friction, which has been found to have a significant effect on the mode propagation [4]. 

The paper is structured as follows. In Section II we present the formalism for studying the oscillations of a two-fluid 
neutron star with mutual friction and vortex lattice elasticity. In Section III we perform a plane wave analysis of the 
oscillation spectrum in the incompressible case and in Section IV we present the more realistic compressible case. As 
we shall see compressibility and chemical coupling between the components can profoundly alter the nature of the 
Tkachenko waves, leading in some cases to much shorter periods of oscillation, close to the rotation period of the star, 
which are not consistent with the observed periodicities in pulsar timing residuals. Finally in Section V we outline 
our conclusions. 



II. TWO FLUID EQUATIONS OF MOTION 

Our starting points will be the multi-fluid formulation of superfiuid hydrodynamics of Andersson and Comer [17] , 
and the Baym-Chandler formalism for including the effects of vortex lattice elasticity in the study of superfiuid 4 He 
|14j . Let us consider a two fluid system of neutrons and protons (which we assume locked to the electrons [TB]) and 
write the Euler equations for the neutrons, in a frame rotating with the star at fixed angular velocity f2 and in the 
absence of external forces and mutual friction. Following [T7] this takes the form: 

(d t + v?V j )(v? + e n wf n ) + 2e ijk Wv* + V,(/i n + cj>) + e n ^ n V lt) ° = (1) 

Where vf is the neutron velocity, wf n — vf — vf, with vf the proton velocity, e n is the entrainment parameter for the 
neutrons [20], An is the chemical potential per unit mass of the neutrons and <j> is the gravitational potential. Note 
that we assume summation of repeated indices and assume the neutron and proton masses equal, m n = m p = m. In 
the above equation we have not yet imposed that the neutrons be superfiuid. To do this we must require that the 
fluid rotates by forming an array of singly quantized vortices, and that averaging over a large number of such vortices 
gives rise to the macroscopic vorticity w ! of the fluid via the relation: 

to 1 = Ku v k l = -e^Vj-Og + e n w p k n ) (2) 

where k % is a unit vector along the direction of the vortex array, K — h/2m n = 1.99 x 10~ 3 cm 2 s _1 is the quantum 
of circulation and n v is the number of vortices threading a unit surface. It is important to remark here that the 
quantization condition on the circulation is a condition on the momentum of the neutron fluid pf = m(vf + e n wf n ), 
and not on its velocity vf (which will not in general be aligned with pf due to the entrainment). From equation |2]) 
one can derive the equation of motion for the circulation 

dtUi + e^e^VW™ = (3) 

and a conservation equation for the vortex number 

a t n v + Vi(n v <) = (4) 

with v l v the macroscopically averaged vortex velocity. It is possible to show ([UJ) that, in order for equations |3| and 
Q to be staisfied it is necessary to add a "Magnus force" term to the right hand side of equation Q, which thus 
takes the form 

(fit + «"V J )« + e n wf n ) + 2 ( , ;i <>•'/•;; + Vi(/t n + <p) + e n w pn W t vJ = Kn v e ijk k j (v* - t£) (5) 

It is clear from equation ^ that in the absence of other forces the vortices will be forced to move with the superfiuid 
neutron condensate. The presence of vortices will, however, also affect the proton fluid, which will experience a drag 
force of the form p n Kn v lZ(vl — v p ), where the exact nature of the process giving rise to the drag is encoded in the 
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dimensionelss parameter 1Z. In a neutron star there are, in fact, a variety of mechanisms that can produce a dissipative 
drag: scattering of electrons off vortex cores is likely to be the dominant process in the core (22j [23] , while in the crust 
the main contribution is due to interactions with the lattice phonons |24| and the excitation of vortex Kelvin waves 
[231 |2"6] . The diverse nature of these processes leads to the drag parameter spanning several orders of magnitude in 
the different regions of a neutron star interior (from as low as 7Z « 10 -10 to 1Z ps 1). We shall thus treat 1Z as a free 
parameter and investigate how its variations affect the modes. 
The Euler equations for the proton fluid take the form: 

(d t + v?V j )(v? - e p wf n ) + 2e ijk Wv k p + V^ P + 0) + £ P <pV^ p = «n v (1 ~ Xp) ^« - vf) (6) 

X p 

where x p = p P /(p n + p p ). p p and e p are now the chemical potential per unit mass and entrainment parameter of 
the protons, such that e p = s n (l — x p )/x p |20) . We also need an equation of motion for the vortex lines which, if we 
assume that they have negligible inertia, takes the form of a force balance between the Magnus force, the drag force 
and the elastic force exerted by the lattice [Ti] : 

p n Kn v e ijk k j («* - «n) + PnKn v TZ(vf - vj) - p n a t = (7) 

where Oi represents the contribution due to lattice elasticity and takes the form: 

2V^(Vi ej ) - (VDeJ (8) 



Pv 

Oi = — 
Pn 



where is the displacement of the vortex line from its equilibrium position, is the gradient perpendicular to 
the direction of the array and p v = p n K 2 n v /16ir is the shear modulus of a triangular vortex lattice [10) . Note that 
the above expression only describes the linear order corrections in the lattice displacements, which are assumed to 
be small. Furthermore we are neglecting the contribution of vortex bending, which would give rise to Kelvin waves 
propagating along the vortex lines. Note that this could be accounted for by including a vortex "tension" term in cr,, 
which we denote af , of the form: 



p n K 2 n v ( b\ d 2 e 



8n \a J dz 2 



where a is the vortex core radius, b is the inter-vortex spacing for a triangular lattice and the z axis is taken along 
the rotation axis of the star. 

The continuity equations for neutrons and protons take the form 

a t p„ + v l (p„<) = (10) 

d t p P + V\p P v?)=0 (11) 

and the gravitational potential obeys the Poisson equation 

V 2 4> = AnG(p n + p p ) (12) 

where G is the gravitational constant. Finally to solve the problem we need to supply an equation of state for the 
system. As we shall examine different cases we delay the discussion of the equation of state to the following sections 
and move on to discussing perturbations of the multi-fluid equations of motion presented above. 

A. Perturbations 

In order to keep the problem tractable we shall consider linear perturbations of a background in which the two 
fluids rotate together with uniform angular velocity Q. For such a background equation ([2]) takes the form: 

ku v = 2il (13) 

and the perturbed Euler equations can be written, in a frame co-rotating with the star, as: 

d t (Svf + e n Swf n ) + 2e ijk tt j 5v* + V^An = -ffiR(6v? - Svf) - a, (14) 

d t (Svf - e p Swf n ) + 2e ijk Sl j 6v* + V t 5~p p = ~ Xp) K(6vJ - Svf) (15) 
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where we have made the Cowling approximation, i.e. neglected the perturbations of the gravitational potential 5(f). 
Note that as a consequence of the extra elastic term in the force balance equation for the vortices (17), the forces on 
the right hand side of the Euler equations are no longer symmetric and the vortex elasticity term only acts on the 
neutron superfluid. The clastic force <jj can be written as 



Crp 



2V^(Vi £i )-(Vi)e, 



(16) 



where we have defined the Tkachenko wave speed c 2 ^ = nfl/Sir and are assuming the vortices to be in equilibrium in 
the background, such that e^ KG = 0. We assume all vortex displacements to be perturbed quantities, and write €i 
in place of Set, unless otherwise specified, and thus consider ai to also be a perturbed quantity. As we are dealing 
with linearised elasticity and the displacement vectors it would be natural to consider Lagrangian perturbations 
of the two-fluid equations of motion, given that in general one would have that Aw* = dt€\. However, given that we 
are working in a rotating frame, and have assumed that the fluids (and thus the vortices) are moving together in the 
background, one has that At;* = Sv^. We can thus continue to work with Eulerian perturbations, which simplifies 
somewhat the problem. 

The equation of force balance for the vortices ^ can be cast in the form: 



SvJ = Svf + 



n 



l+n i e m- "~pn 2 ft(l+ft 2 p 1 l + TZ 2 ^" 1 "' 1 

and the perturbed continuity equations, in the absence of reactions, take the form 

d t S Pn + V i ( Pn Svf) = 
d t 5p p + V i ( Pn 5vf) = 



n 



f.K J 0W 



;(5w 



8w l - 



K 2 2n(l+K 2 



Following [27] we can combine equations ( 15 ) to obtain an Euler equation for the "total" velocity Vi 



d t Svi + -ViSp - 
P 



-Vjp + 2e ijk Q, :l 5v =-(l-x p )ai 



(17) 

(18) 
(19) 

-x p )vf+x p vf: 
(20) 



and one for wf n : 



(1 - s)dtwf n + V t 5p = -2Q.B' e m k j Swl n + 2Q,Be ijh k? e klm ki8w^ 



where we have defined the total pressure, such that V^P = Pn^iP-n + Pp VjjEt p , the entrainment parameter e 
6/3 = Sp p — Sp n and the mutual friction parameters B 
continuity equations ( 19 1 can be cast in the form: 



(21) 



l-K 2 /[x p (l+K 2 )] and 6 = K/[x p (l + K 2 )]. The perturbed 



d t 5p + V t { P 5v l ) =0 
d t Sx p + -Vj[x p (l - x p )pSi 



Sv-'VjXp = 



(22) 
(23) 



As we shall see in the following, this formulation can be advantageous when discussing the compressible problem. 



III. THE INCOMPRESSIBLE CASE 



In order to make contact with previous results, let us consider first of all the case of incompressible fluids, such that 
Sp = and the continuity equations reduce to 

X7*5vf = V*S< = (24) 

We consider plane waves, such that a perturbed quantity 5/,(x, t) takes the form 5/j(x, t) — fi exp(ikiX l — icut), with 
fi a constant amplitude. Without loss of generality we choose our coordinate system such that the z axis points along 
the rotation axis and such that the wave vector k % lies in the x — z plane, i.e. k = (fcsin#, 0, kcos9). The equations 
of motion can thus be written as: 

— iu)vf(l — e n ) — icje n vf + 2eijk^v k + ikifi n = 2WZ(iujei + vf) — <jj (25) 
-iu)vf(l - e p ) - iuje p vf + 2e ljk n 1 v k + ik % p p = -20 ^ ~ Xp ' H(iuei + vf) (26) 

kj4 = kj4 = (28) 
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where we have defined B = 11/(1 + TZ 2 ), di = ai/20, and to simplify notation we have defined /i x as the amplitude of 
5jl x , with x — n, p. Finally the vortex elasticity contribution takes the form: 

er = qe with q = (-(c T ksm8) 2 , (c T ksin8) 2 , 0) (29) 

In order to obtain the dispersion relation for the modes of the system we thus need to solve the characteristic equation 
det|.Ky|=0, where follows from equations ( 25]p8 l and is given in equation (Al). 



In the undamped case, neglecting the effect ofentrainment (e n = e p = 0) one obtains, as expected, two families of 
modes, the inertial modes 

u 2 = 4fl 2 (cos8) 2 (30) 

and the Tkachenko waves 

lj 2 = 4ft 2 (cos(9) 2 +4fc 2 (sin(9) 4 - -^-(sinfl) 4 « 4tt 2 | cos6»| 2 + c 2 T k 2 {sm9f (31) 

where we are assuming that c^k 2 << fl 2 . This will always be the case if we consider typical pulsar spin rates 
from a few Hz to a few hundred Hz and long wavelength oscillations across the whole superfluid region, such that 
k s» 10~ 5 ~ 10~ 6 cm -1 . For propagation perpendicular to the rotation axis (cosO = 0) one then obtains the well 
known Tkachenko wave dispersion relation 

uj = ±c T k (32) 

A. The effect of entrainment 

Let us now still consider undamped propagation of the modes, but include the effect of entrainment. Clearly 
introducing coupling between the two fluids profoundly alters the nature of the modes and leads, in the c^fc 2 << il 2 
limit, to two families of mixed inertial-Tkachenko waves: 

w 2 4ft 2 (cos6>) 2 + (1 -x p )c 2 T k 2 (sui8) i (33) 
lu 2 » 4» 2 (cos0) 2 f— gg— V +c 2 T k 2 j (sing) 4 (34) 

In the limit x p — > I and e n — > one has again e p = e n = from the relatio n e p = e n (l — x p )/x p , the two fluids 



decouple and we have the two separate families of modes of equations (31) and (32) 



B. Mutual friction 

We now consider the dissipative terms due to mutual friction, i.e. to the drag parameter TZ. In order to keep the 
results tractable we take e n = e p = 0. It is still however impractical to consider the whole solution, so let us first 
of all consider modes propagating along the z axis (8 = 0). In this case one has two families of inertial modes, one 
which is undamped with dispersion relation 

ui = ±2n (35) 

and one which is affected by mutual friction: 

uj = ±2flB - i2flB (36) 



where we remind the reader that B = 1 - K 2 /[x p (l + TZ 2 )} and B = H/[x p (l + TZ 2 )]. The results in (|35j and (|36|) 
agree well with those of [27], in which the authors show that there is one class of inertial modes that corresponds 
to the fluids co-moving and is undamped (in the absence of chemical coupling) and another class of counter-moving 
modes that is rapidly damped by mutual friction. 

Let us now examine the case of Tkachenko waves propagating perpendicular to the rotation axis (cos 8 = 0). In 
figure [l] we plot the frequency of the modes for k = 10~ 6 and f s tar = 10 H z i as a function of TZ in the weak drag 
regime. For large values of the proton fraction x p we recover the solution of in which the real part of the frequency 
vanishes and the damping becomes large for values of the drag parameter such that the mutual friction damping 
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FIG. 1: We plot the real part of the modes and the modulus of the imaginary part (dotted lines), for x p =0.96 and x p =0.1. 
We take k = 10" 6 . For x p = 0.96 we recover the results of [4] : the real part of the mode vanishes when the mutual friction 
damping timescale is close to the mode period, and there is an extra purely imaginary root. For the more realistic, but still 
large, value of x p —0.1 we see that, on the other hand, the mode is always oscillatory, although the imaginary part is larger 
when the damping timescale and mode period are similar. For higher values of the drag the frequency of the mode is reduced 
to « 25% of the original value. 



timescale r m ~ 1/2QTZ is approximately equal to the the Tkachenko wave period Pt = 2tt /lot with ujt = ky/nQ/n. 
In fact close to this value the mode has an extra purely imaginary root, as found by [3]. For more realistic values 
of the proton fraction we see, however, that the pathological behaviour disappears and even though the damping is 
stronger when the mutual friction timescale is close to the period of the modes, the real part does not vanish, the 
mode is always oscillatory, and there are always only two purely imaginary roots. 



C. Perfect pinning 



Up to now we have assumed that the vortex lines are free to move and experience a drag force as they do so. 
However it is commonly believed that vortex lines can interact strongly with lattice points in the neutron star crust 
and 'pin' to them, in such a way that they are forced to move with the charged components of the star [28H32] . The 
nature of such a pinning force is well beyond the scope of this paper, but to study the propagation of Tkachenko 
waves in this scenario it is sufficient to consider an unspecified force / pin acting on the vortices such that they are 
forced to move with the proton fluid, i.e. such that 

Svl = 6vl (37) 

In this case the equation of force balance for vortices takes the form: 

/'::K'M, ;i ./r< (Sv k p - 8v*) - Pn a t + r in = (38) 

where there is no drag force acting, as the vortex lines flow with the protons. If we now consider the neutron and 
proton fluid there will be a reaction force —f pin acting on a the protons and the Magnus force acting on the neutrons. 
Making use of equation (38) we can thus cast the Euler equations in the form: 

d t (5vf + e n 6wf a ) + 2e ijk W6v k + V^fc = -2Sle ijk k j (6v* - Sv k ) (39) 

d t (5v? - e p Swf n ) + 2e ljk W6v k + V i( 5/i p = 2ft (1 ~ Xp h ijk K j (5v* - Sv k ) - ^ ~ Xp) a { (40) 



which lead to: 



-iuvf(l - e n ) - iuJE n barvf + 2e ijk n j v^ + ik^ n = -2Q,e ijk k j {v^ - v k ) (41) 
-iw»?(l - e p ) - iwe p «? + 2e ijk Wvl + ikifiv = 2n ( - 1 ~ Xp K ijk k :i (v^ - v k n ) - (1 ~ Xp) a, (42) 
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The characteristic equation given by the equations in ( |42|), together with the condition in (37) can be obtained by 
calculating the determinant of the matrix Kij given in ( |A12||A18 1. This leads to two families of modes: 



uj 2 = =4tf(cosd) 2 ( L^A +4fc 2 ( sm fl) 4 (1 X *f (43) 



,2 



to 



= = 4fi 2 (cos6») 2 + c|fc 2 (sin6') 4 (l - x p ) (44) 



These are once again mixed inertial-Tkachenko waves, but we see that in the limit x p — > 1 the Tkachenko waves 
disappear and we are left with only one family of inertial modes. This resembles the situation in superfluid 4 He, in 
which even a small amount of pinning swamps the contribution due to lattice elasticity and transforms the Tkachenko 
waves into inertial waves |33j . 

IV. COMPRESSIBLE NEUTRON STAR MATTER 

It is well known from the study of superfluid 4 He that compressibility can have a drastic effect on the mode structure 
[36] . Including compressibility in the equations of motion for the superfluid leads to a dispersion relation of the form 
[33] : 

„2 .2 1.4 

-'^^y < 45 > 

where c s is the sound speed. In the long wavelength limit (fc << f2/c s ) the nature of the mode is thus profoundly 
altered and the dispersion relation is no longer linear in fc, but rather parabolic, leading to the so-called "soft" 
Tkachenko wave frequency: 

In the study of 4 Hc the long wavelength limit is, however, mainly of theoretical interest, as one would need containers 
of several hundreds of meters in diameter to explore it experimentally. The situation is very different for BECs as, 
in contrast with a strongly interacting Bose liquid such as He, they are weakly interacting Bose gases with low 
sound speeds for which the effect of compressibility is important at high rotation rates. For BECs the compressible 
Tkachenko wave spectrum has thus been studied both theoretically [T21II31135 and experimentally [33] , 

Let us now consider a realistic neutron star . The situation is clearly quite complex as not only can we be in the 
long wavelength limit (f2 « c s fc) for the more rapidly rotating pulsars, but one also has to account for multi-fluid 
effects and chemical coupling between the different constituents via the equation of state. One cannot, in general, 
assume incompressibility for the proton and neutron fluids and it is clearly of great interest to adapt our formalism 
to include the effects of compressibility and chemical coupling. To study this problem it is now advantageous to write 



the perturbation equations in the form of equations (20)- (23). In the plane wave approximation the Euler equations 
take the form (in the Cowling approximation): 

—iuvi + i—p - -ViP + 2e ijk fl j v k = -(1 - x p )<7 t (47) 



P 



K 



—iuj(l — e)wi + ih/3 + 2eijk^l J w — — 2Q — (iwe; + Vi + (1 — x p )wi) + er,; (48) 

x p 

and the continuity equations can be written as: 

-iujp + ipkjV j + v 3 Vjp = (49) 
— iuix p + ix p (l — x p )kjW 3 + w ] Vj \px p (l — x p )] — 0, (50) 



while the equation of force balance still takes the form in (27). As we are now considering compressible matter, we 



will also need an equation of state for the perturbations. Choosing to work with the density (p) and proton fraction 
(x p ) perturbations, one can write: 

*Vf#-U (51) 

(52) 



dp J \ dx p 
<9/3\ _ / d/3 



dp J \dx 
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52) should be derived from a fully 



Ideally the partial derivatives of the thermodynamical variables in equation 
consistent multi-parameter equation of state, which should also allow us to calculate the entrainment parameters and 
the superfluid gaps for neutrons and protons. However, not only is such an equation of state not currently available, 
but its use would also be beyond the scope of our simplified plane wave analysis. In order to keep the problem tractable 
we shall assume that our background model is described by an n — 1 polytrope and use for the perturbations two 
simplified analytic equations of state, that are essentially extensions of a single-fluid n = 1 polytrope. 

First of all we shall consider the equation of state of [15], which we refer to as model A. In this case we have: 

PA (53) 
c 2 

= d (54) 

where c s is the sound speed of the background We shall then consider a second model, in order to understand the 
importance of the chemical coupling on the mode structure. This model, which we refer to as model B, is essentially 
a re-parametrisation of the model II equation of state in [37J (also used in [SH] where it is denoted as model BO) and 
takes the form: 









dp) 




\dx p ) 




= ^-, 




dp) 


pxp 1 


\dx p 



dp\ _ ( dp \ pc\ 



dp) ° sl \dx p ) a Xp ^ 

dp) a pxp \dxp) x\ 

We can thus study the behaviour of the solutions to are problem as we vary the parameters a and 7. Clearly model 
A corresponds to the case a = 1, 7 = 1. 

The sound speed in the background, for ann = l polytrope, takes the form 

c s = 2Kp (57) 

where K = 2GR 2 /n depends only on the stellar radius. However in our plane wave approximation we shall assume 
that the background quantities vary over a length-scale greater than that of the oscillations, and thus take them to be 
constant and neglect their gradients. This approximation is not necessarily justified, as in the crust the density and 
pressure vary by several orders of magnitude over a length-scale of approximately 1 km, which is comparable with 
the longest wavelengths we consider for our Tkachenko waves. It is, however, a reasonable approximation for shorter 
wavelengths and in the neutron star core. The sound speed will thus be a constant in our formulation and specifically 
we take c s = 10 9 cm s _1 . We also take the proton fraction, which in a rigorous description should also be derived 
from the equation of state, as a constant and will study the effect that varying it can have on the modes. Needless to 
say, future work should aim to relax this approximation and consider a fully stratified neutron star. 

Finally let us remark that for simplicity we take e = in the following discussion. We have experimented with 
varying the parameter e between -0.8 and 0.8, but it is found to have very little effect on the dispersion relation. 



A. Undamped propagation 

1. Model A 

Let us consider, first of all, the undamped propagation of waves in a neutron star. We thus take 1Z — and, to 
keep the problem tractable, e — 0, and solve the characteristic equation obtained form the determinant of the matrix 



Kij given in (Bl I. As a first step we focus on model A. The results are two families of modes, such that: 



u? = 40 2 - 4fc y 4 (58) 
u 2 = ±\ (An 2 + k 2 c 2 s ) + ^-J{k 2 c 2 s + in 2 x p ) 2 - Xp(4flkc s cos0) 2 (59) 



As we can see we obtain sound waves and rotationally corrected Tkachenko waves, for which the main contribution 
to the frequency is now given by the stellar rotation frequency. This result is somewhat surprising, as in this limit 
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FIG. 2: We plot the frequencies of the modes we obtain, normalised to the classical Tkachenko wave frequency, for two different 
rotation rates of the star and for a varying parameter 7, whilst keeping a = 0. We take x p — 0.05. We can see that we have two 
families of high frequency sound waves and then the Tkachenko waves, the frequency of which is shifted from the classical value 
at higher rotation rates. This is expected as for higher rotation rates the effects of compressibility become more important. 
Furthermore it is clear from the graph that, although varying 7 has a small effect on the frequency of the sound waves, it has 
no effect on the frequency of the Tkachenko waves. In these plots we have set e — 0, k — 10 -6 and taken 7 = 1.5. 



the classical Tkachenko waves no longer exist and the vortex elasticity simply provides a small correction to what 
is, in essence, the frequency of an inertial wave. In this case, for a pulsar rotating at ~ 1Hz, the frequency of the 
Tkachenko waves would be much too high to explain the observed periodicities of 100 ~ 1000 days observed in the 
timing residuals. Such a drastic modification in the dispersion relation clearly needs to be investigated in more detail. 
We have, after all, used a simplified version of the equation of state, so let us turn our attention to model B in order 
to understand how varying the parameters (and thus the coupling between the components) can affect the mode 
structure and whether there is a reasonable set of parameters for which one can still obtain the usual Tkachenko 
waves. 



2. Model B 



In the case of model B the extra parameters make it necessary to solve the characteristic equation numerically, but 
they also allow us extra freedom to explore different regimes for the equation of state. First of all let us eliminate 
all chemical coupling by setting a = 0. In figure [2] we plot the mode frequency for varying 7 at different rotation 
rates. We see that in this case there is a mode with the frequency of the classical Tkachenko mode and furthermore 
varying 7 has very little effect on its frequency. If we now keep 7 fixed and vary the rotation rate of the star we see 
in figure [3] that for low rotation rates the frequency of the Tkachenko mode is the classical one. For higher rotation 
rates the classical frequency increases but the frequency of the Tkachenko mode tends to that of the so called "soft" 
Tkachenko mode CyC s fc 2 /20. The result is thus that, for longer wave-length of order the stellar radius, there is 
always a Tkachenko mode with periods consistent with the 100 — 1000 day variability typical of pulsar timing noise. 
This picture is clearly very different from that of the previous section, in which the Tkachenko waves had essentially 
disappeared, so let us investigate how re-instating the chemical coupling and varying the parameter a can affect the 
mode structure. 

In figure [4] we plot the mode frequency for a stellar rotation frequency of 10 Hz varying a whilst keeping 7 fixed. 
The result is now much more intriguing as one still has one family of sound waves, but there is then an avoided 
crossing between the second family of sound waves and the Tkachenko waves, with the frequency of the sound wave 
becoming that of a classical Tkachenko wave as we vary a. The real part of frequency of the Tkachenko wave, on the 
other hand, vanishes after the avoided crossing and one obtains two purely imaginary roots. 

Summarising there is a vast portion of parameter space in which one has a mode close to the frequency of a classical 
Tkachenko wave but there exists a small region (which thus includes model A for which a — 7 = 1) where the avoided 
crossing occurs, in which the Tkachenko mode is not oscillatory in nature and the frequency of the sound waves is too 
high to account for the slow variability of pulsar timing residuals. If we now increase the stellar rotation rate we can 
see from figure ^ that the region in which the Tkachenko mode is not oscillatory becomes even larger, allowing for 
vast portions of parameter space in which there is no mode close to the Tkachenko frequency. However we have not 
yet considered the impact of mutual friction damping, which could potentially limit even more the range of parameters 
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FIG. 3: We plot the frequencies of the Tkachenko waves we obtain numerically and compare them to the classical Tkachenko 
wave frequency crk and to the "soft" Tkachenko wave frequency CTC a k 2 /2Q. As we can see the frequency tends to that of the 
"soft" mode for higher rotation frequencies (in the millisecond range, which is that of the fastest known pulsars) and is slightly 
modified by multi-fluid effects for high values of x p . In these plots we have set e — and taken k = 10~ 6 . 




FIG. 4: For a stellar rotation rate of 10 Hz, we plot the frequency of the 
frequency, for varying values of the parameter a. We see that there is still 
is now an avoided crossing between the second family of sound waves (p2) 
of the sound wave becoming that of a classical Tkachenko wave. There 
one has a mode close to the frequency of a classical Tkachenko wave, but 
occurs, in which the Tkachenko mode is not oscillatory in nature and the 
for the slow variability of pulsar timing residuals. Once again we have set 



modes, normalised to the classical Tkachenko mode 
a family of sound waves (indicated as pi), but there 
and the Tkachenko waves (Tk), with the frequency 
is thus a vast portion of parameter space in which 
there exists a small region where the mode crossing 
frequency of the sound waves is too high to account 
e = and taken k = 10 -6 . 



for which one has long lived Tkachenko oscillations. Let us thus move on to consider the full problem. 



B. Mutual Friction 



The inclusion of mutual friction makes the problem considerably more complicated, and once again the characteristic 
equation must be solved numerically. The picture that emerges is however not very different from that of the previous 
section. One still finds a high frequency family of sound waves, largely unaffected by mutual friction, then a second 
family of sound waves at lower frequency and a family of Tkachenko waves. In addition one has two purely imaginary 
roots to the characteristic equation. Let us focus on the Tkachenko waves and on the lower frequency sound waves. 
As we can see in figure ([6| the behaviour of the modes depends strongly on the chemical coupling. For (3 > a and 
slow rotation of the star one has Tkachenko waves close to the classical frequency, while for j3 < a one has an avoided 
crossing between the second family of sound waves and the Tkachenko waves. In both cases the waves oscillating close 
to Tkachenko frequency are strongly damped by mutual friction in a narrow range of the parameter 1Z for large values 
of a; p , exactly as in the incompressible case. For more realistic values of x p we find that, as expected, the damping is 
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FIG. 5: We plot the frequencies of the modes we obtain (sound waves pi and p2 and Tkachenko waves Tk), normalised to the 
classical Tkachenko wave frequency, for two different rotation rates of the star and for a varying parameter a, whilst keeping 
7 = 0.8. We take x p — 0.05. We can see that as the rotation rate increases not only does the Tkachenko wave frequency 
decrease as expected, but there is also a vast region of parameter space in which the Tkachenko mode dissapears. In these plots 
again we have set e = and taken k = 10 -6 . 



negligible. 

If we now increase the rotation rate in the case (3 > a one finds that the frequency of the Tkachenko waves 
approaches that of the "soft" mode and mutual friction only weakly damps the mode for realistic values of x p . The 
picture is considerably different if we take /3 < a, as can be seen from figure Q. For a rotation rate of 60 Hz one 
has a sound wave close to the classical Tkachenko wave frequency and a highly damped "soft" Tkachenko wave for 
low values of 7Z. However, if we increase the rotation rate to 100 Hz, the frequency of the "soft" mode becomes 
purely imaginary and the sound waves return to being high frequency modes, weakly damped by mutual friction. 
The picture that emerges is thus that, while for low rotation rates of order of a few Hz one has long lived Tkachenko 
oscillations for a vast range of plausible parameters (except for the particular case of a = f3 — 1 as in model A), for 
higher rotation rates, approaching 100 Hz, the situation is radically different and for several choices of parameters 
there are no Tkachenko modes at all. 



C. Perfect pinning 

Finally we examine, as in section III C, the case of perfect pinning. The Euler equations can be written as: 

-iujVi + i-pP- ~^iP + 2e ijk n : >v k = -(1 - x p )(7i (60) 
- i)uw t + ikS - 2 ( - 1 ~ Xp K l]k Ww k = - ^~ Xp K i (61) 

Xp Xp 

Together with the pinning condition —iujti = vf = Vi + (1 — x p )wi. In the limit of no entrainment (e — 0) for purely 
transverse propagation (cos6> = 0), one finds that, as in the incompressible case, the spectrum depends heavily on the 
value of the proton fraction x p . For slow rotation rates (of the order « 10 Hz) one always has a Tkachenko mode for 
low values of the proton fraction (less than x p « 0.3). However the situation changes for higher rotation rates. For a 
stellar rotation rate of 100 Hz, one can see from figure ^ that for the, possibly more realistic, case of small values of 
Xp one has a Tkachenko mode only for 7 > a. However for larger values of a; p the opposite is true and the Tkachenko 
mode only exists in the limit 7 < a. If 7 = a = 1 (model A) one has, as expected, no Tkachenko waves for any value 
of the proton fraction. 



V. CONCLUSIONS 



In this paper we present a formalism for the inclusion of vortex lattice elasticity in the multi-fluid hydrodynamics 
of superfluid neutron stars. As a first step we consider incompressible neutron and proton matter, in order to make 
contact with the well known results for superfluid 4 He. We obtain the standard dispersion relation for Tkachenko 
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waves and find that, in the limit in which p n ~ Pp, mutual friction completely damps out the oscillations in less than 
a period when the damping timescale due to the drag is close to the mode period (as found by However in a 
neutron star one will, in general, have that p p << p n , and in this case we always find weakly damped oscillatory 
solutions. It should be stressed that a more realistic model should include other sources of damping, such as shear and 
bulk viscosity, and the calculation should be performed for global modes in spherical symmetry, without assuming 
constant background quantities. 

The main focus of this work is on the effect of compressibility and chemical coupling on the mode spectrum. We 
find that for slow rotation rates (a few Hz) one can, in general, obtain solutions that correspond to low frequency 
Tkachenko waves and could explain the timing noise in older pulsars. However, for particular choices of the EOS 
parameters, such as model A, there are no Tkachenko waves, but only modified inertial waves and sound waves, 
that oscillate at frequencies that are to high to have any connection with the timing noise. Finally the situation is 
considerably more complicated for more rapidly roatating NSs (above 100 Hz) for which there is now a large portion 
of parameter space for which there are no propagating Tkachenko waves. It is thus clearly imperative to obtain more 
stringent constraints on the EOS from nuclear physics, in order to understand if the regions of parameter space in 
which one has no Tkachenko waves are of physical significance or not. In the light of these uncertainties, it is still very 
much an open question whether or not the period of Tkachenko waves in a realistic neutron star could explain the 
observed periodicity of « 256 days in the timing residuals of PSR B1828-11 (which is rotating at v — 2.469 Hz) and 
the timing noise in other pulsars, or if it is likely to power low frequency precessional motion of the star. Our results, 
however, indicate that for a large range of physical parameters long period Tkachenko waves can in fact propagate in 
a NS interior and are likely to play a role in the dynamics of the system. 

Finally let us remark that we have considered perturbations of a co-moving background. While this is not a bad 
approximation in many situations, it is possible that if the vortices are pinned to the crust a significant lag could 
build up between the charged component and the superfluid neutrons. This can lead to a series of short wave-length 
instabilities ([39], [40]) that are likely to have an impact on pulsar glitches and on the stellar response to external 
torques, such as those experienced by neutron stars in accreting systems. We plan to relax the assumption of a 
co-moving background and explore the consequences on these physical scenarios in future work. 
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Appendix A: Characteristic equation: the incompressible case 
For the incompressible case the equations of motion ( 25|28 ) can be cast in the form: 
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and we recall that the displacement e J only has components in the plane perpendicular to the vortices, i.e. the x — y 
plane in our formulation. The components of the matrix Kij take the explicit form: 
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where we recall the definition B = 1Z/(1 + 1Z 2 ). The characteristic equation then follows from the determinant of the 
matrix Ka. 



1. perfect pinning 

In the case of perfect pinning the elements of the matrix are modified in the following way: 
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1. perfect pinning 

In the case of perfect pinning, for a compressible model, the elements of the matrix Kij are modified in the following 
way: 
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